{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Need to figure out how to integrate the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "TARGET_RADIUS = 1.5 * 25.4 # mm\n",
    "TARGET_TILT = 15 * np.pi / 180 # rad\n",
    "SUBSTRATE_RADIUS = 10 # mm\n",
    "TARGET_SPIN_RATE = 1 # per second\n",
    "\n",
    "CHAMBER_VERTICAL = 136\n",
    "CHAMBER_LATERAL = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cart_to_polar(x, y):\n",
    "\n",
    "    r = np.sqrt(x**2 + y**2)\n",
    "    phi = np.arctan(y/x)\n",
    "    return r, phi\n",
    "\n",
    "def cylinder_sag(y, RoC):\n",
    "    return -np.sqrt(RoC**2 - y**2) + RoC\n",
    "\n",
    "def cylinder_surface_normal(y, RoC):\n",
    "    denomenator = np.sqrt(RoC**2 - y**2)\n",
    "    xcomp = np.zeros_like(y)\n",
    "    ycomp = y / denomenator\n",
    "    zcomp = np.ones_like(y)\n",
    "\n",
    "    vec = np.array([xcomp, ycomp, zcomp])\n",
    "\n",
    "    if vec.ndim > 1:\n",
    "        for _ in range(vec.ndim - 1):\n",
    "            vec = np.moveaxis(vec,-1,0)\n",
    "\n",
    "    return vec\n",
    "\n",
    "def erosion_profile(a, max_radius):\n",
    "    return np.sin(a * np.pi / max_radius)\n",
    "\n",
    "def target_surface_normal(omega):\n",
    "\n",
    "    xcomp = -np.sin(omega)\n",
    "    ycomp = 0.\n",
    "    zcomp = -np.cos(omega)\n",
    "\n",
    "    return np.array([xcomp, ycomp, zcomp])\n",
    "\n",
    "def distance_to_target(x, y, a, phi, h, l, omega, zfunc, znorm, RoC, spin_rate, t, spin_offset=0):\n",
    "    \"\"\"compute distance from point on target to a point on the substrate\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : float\n",
    "        cartesian point on the substrate\n",
    "    y : float\n",
    "        cartesian point on the substrate\n",
    "    a : float\n",
    "        radial point on the target\n",
    "    phi : float\n",
    "        azimuthal point on the target\n",
    "    h : float\n",
    "        vertical distance from target to substrate\n",
    "    l : float\n",
    "        lateral distance from target to substrate\n",
    "    omega : float\n",
    "        tilt of the target w.r.t. the horizontal axis\n",
    "    zfunc : callable\n",
    "        sag function of the substrate, i.e. sag = zfunc(y, RoC)\n",
    "    znorm : callable\n",
    "        surface normal function of the substrate, i.e. normal = znorm(y, RoC)\n",
    "    RoC : float\n",
    "        Radius of curvature of the substrate\n",
    "    spin_rate : float\n",
    "        angular frequency of the substrate rotation\n",
    "    t : float\n",
    "        time in seconds\n",
    "    spin_offset : float\n",
    "        substrate angular offset from alignment to the target axis of lateral displacement\n",
    "    \"\"\"\n",
    "\n",
    "    r = np.sqrt(x**2 + y**2)\n",
    "    cosphi = np.cos(phi)\n",
    "    sinphi = np.cos(phi)\n",
    "    cosomega = np.cos(omega)\n",
    "    sinomega = np.sin(omega)\n",
    "\n",
    "    z = zfunc(y, RoC)\n",
    "    ns = znorm(y, RoC)\n",
    "    ns /= np.sqrt(np.sum(ns * ns, axis=-1))[..., np.newaxis]\n",
    "\n",
    "    nt = target_surface_normal(omega)\n",
    "    nt /= np.sqrt(np.sum(nt * nt))\n",
    "\n",
    "    xcomp = r * np.cos(spin_rate*t + spin_offset) - l\n",
    "    xcomp -= a * cosphi * cosomega\n",
    "\n",
    "    ycomp = r * np.sin(spin_rate*t + spin_offset) - a * sinphi\n",
    "\n",
    "    zcomp = a * cosphi * cosomega \n",
    "    zcomp -= h - z\n",
    "\n",
    "    vec = np.array([xcomp, ycomp, zcomp])\n",
    "\n",
    "    if vec.ndim > 1:\n",
    "        for _ in range(vec.ndim - 1):\n",
    "            vec = np.moveaxis(vec,-1,0)\n",
    "\n",
    "    return vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Target geometry\n",
    "u = np.linspace(-TARGET_RADIUS, TARGET_RADIUS)\n",
    "u, v = np.meshgrid(u, u)\n",
    "a, phi = cart_to_polar(u, v)\n",
    "\n",
    "# Substrate geometry\n",
    "x = np.linspace(-SUBSTRATE_RADIUS, SUBSTRATE_RADIUS)\n",
    "x, y = np.meshgrid(x, x)\n",
    "RoC = 1e20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = np.linspace(0,100,10)\n",
    "mass_collision = 0\n",
    "\n",
    "for t in times:\n",
    "\n",
    "    # get erosion profile\n",
    "    eps = erosion_profile(a, TARGET_RADIUS)\n",
    "\n",
    "    # get target distance\n",
    "    R = distance_to_target(x, y, a, phi,\n",
    "                           h=CHAMBER_VERTICAL,\n",
    "                           l=CHAMBER_LATERAL,\n",
    "                           omega=TARGET_TILT,\n",
    "                           zfunc=cylinder_sag,\n",
    "                           znorm=cylinder_surface_normal,\n",
    "                           RoC=RoC,\n",
    "                           spin_rate=TARGET_SPIN_RATE,\n",
    "                           t=t)\n",
    "    \n",
    "    norm = np.sum(R * R, axis=-1)\n",
    "    eps /= norm\n",
    "\n",
    "    # get distance direction\n",
    "    Rnorm = R / np.sqrt(norm)[..., np.newaxis]\n",
    "\n",
    "    # get surface normal directions\n",
    "    ns = cylinder_surface_normal(y, RoC)\n",
    "    ns /= np.sqrt(np.sum(ns * ns, axis=-1))[..., np.newaxis]\n",
    "\n",
    "    # get target normal directions\n",
    "    nt = target_surface_normal(TARGET_TILT)\n",
    "    nt /= np.sqrt(np.sum(nt * nt))\n",
    "\n",
    "    target_dot_distance = np.sum(nt * Rnorm, axis=-1)\n",
    "    surface_dot_distance = -1 * np.sum(ns * Rnorm, axis=-1)\n",
    "\n",
    "    mass_collision += eps * target_dot_distance * surface_dot_distance * a\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAGfCAYAAAApoGrxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABNnElEQVR4nO3df5AU5Z0/8Hf3/Nyfo4jMwFciqyH4g9IiGBEqCpeEjZofeiQVDbk9vSQeBn8R6soEuTuX3IUVYlEktYiFMeqVIfiH4WJ9y+yxVxdXU4ABZCMxhKrUEdxvZFzRZXf258x0P98/ViauO8/n2e3Zgen0+1U1RTHPPL09PT29zz7d7/5YSikFIiIiChT7XK8AERERnX0cABAREQUQBwBEREQBxAEAERFRAHEAQEREFEAcABAREQUQBwBEREQBxAEAERFRAHEAQEREFEAcABAREQVQuFwLfuyxx/CDH/wAJ0+exJVXXomtW7fi+uuvN/ZzXRdvvfUW6urqYFlWuVaPiIjKRCmFTCaDWbNmwbbL93fm8PAwstlsycuJRqOIx+NTsEY+o8pg165dKhKJqCeeeEL9/ve/Vw888ICqqalRJ06cMPbt6upSAPjggw8++PD5o6urqxy/YpRSSg0NDanUjNCUrGcqlVJDQ0NlW9dKZSk19cWAFi1ahI9//OPYvn174bnLL78ct956K1paWsS+vb29OO+88/BJ3IwwIuPa7XhM7G/V1Ojb6vRtAODWVmnbnFr9z83Vjl/PD8rW6UfA+Wp5dJyr1c+C5Kr1/RzDYNap0n/sTtwV+6qYvt2KO9q2UFhebjia17ZFDH2jYf3PDdtSm7xc29JvJ0toM/WVuEqe+VJCu6lv3hX2RTekbcvm9W0AkMsLy83KE42O0FcN63+uNSJ/d0LD+vbQkLydQsP6tsig0NYvf+bhQf3+Fs3I+2KkP6dtC/WPaNvs/iFxuSozoG8b0LcBgDtc/OfmkcOv8SJOnz6NRCIhLsOrvr4+JBIJHD90MeqFY6xxORkXDQtPoLe3F/X19VO4hpVvyk8BZLNZHDp0CN/97nfHPN/Y2Ii9e/eOe/3IyAhGRv6yE2UymfdXLIKwVWQAYEXFn2/Z+nbLlgcPbkjfboX1v1FVRB4AuBHhABeVd1w3KhzspbdjGACouP5AparKMwCwI/JyQ8IAICT8gje1h0IcABT6CgMA5eh/2YYMAwBXaHfD8mFG5YR1soQBgGFq2RYucQoZtlNI+OhC+t0UoahhAJDT729hw/cjHBY+H+HjsQ37uLL1Awtl6dsAwLU0y35/M5yN07j1dXZJA4Agm/KtdurUKTiOg2QyOeb5ZDKJdDo97vUtLS1IJBKFx+zZs6d6lYiI6K+Uo9ySH0FVtmHTh0d+Sqmio8F169aht7e38Ojq6irXKhER0V8ZF6rkR1BN+SmA6dOnIxQKjftrv7u7e9ysAADEYjHEYvLUPBERUTEuXJTyN3xpvf1tygcA0WgUCxcuRHt7O/72b/+28Hx7eztuueWWCS/Hrqkuer7fqpUv5INwoZ90kR8A5OuEC/3q9ZtqpN5wgZR0IZ/QBgA54e3ma4QL+YSL/ADArdKfF7eENgCIxvQnQWMx/TnDeEQ4eQogHta3V4Xlc5ExoW/Y8n4NQFh3jhPez/GXSjrPn1fypJ54EaDS78cjeflQMZTXXwczbOg7nBO+W3H9cnMj8nLzMeG6BKENAEIx4TqLiLe20Xb99nfDhus3hFWO2Pq+YcO1ElZI325JFxcAsDXttsoC8vWDVAHKch+AtWvXoqmpCddccw0WL16MHTt24M0338Tdd99djh9HREQB5SgFp4QwWyl9/a4sA4DbbrsN7777Lr73ve/h5MmTmD9/Pl588UVcfPHF5fhxREQUUKWex+c1AGWwevVqrF69ulyLJyIiohKUbQBARERUbi4UHM4AeMIBABER+RZPAXjH2ycREREFEGcAiIjIt5gC8K5iBwB2bS3sIvf1V6aCPvX6rL+U8weAbEK/ObK1+smSbJ0hyy+05w23NchJWf8aIctfLWf5I3F9rr5KaAOA2ri+8Eh1RN+3JiyX7awW2qtChvsACDdol+4DEBHaADnrb7oPgOm+/OVYruln5oQwuXgfAMdwHwBHn9cfzMv1OwaE9kHhPgD9w/L3eSiq75uPyfU7clGptoG+TRmy/FLWX6jFNLpsW1on/XJVSF6nsNBuG+7lr2u33bN3HwD3/Ucp/YOKpwCIiIgCqGJnAIiIiEycElMApfT1Ow4AiIjItxw1+iilf1BxAEBERL7FawC84zUAREREAcQZACIi8i0XFhx4S9yc6R9UFTsAUHU1UKHxMR+3Li72y9fro0FSzA8ARur0MRsp6perFxeLXK3+JJNU0hcAXCHqF6oRyudW66N6gBzlq4/Kfeuiw/q2iL5vXVjfD5CjfnHbEAO09NsiYuvbQoYLgGyhHHAlcg3lgKUDZc4VyvIqQ0lfVx+rkyKCAJDJ67/TmZz++5yJyMeCvpi+b7/QBgBDEX27E9bvM0qICAKlxfWUUPJXajP+fhPaTb8gdDFA5YwAbxs6TxFXjT5K6R9UPAVAREQUQBU7A0BERGTilHgKoJS+fscBABER+RYHAN7xFAAREVEAcQaAiIh8y1WW57obZ/oHFQcARETkWzwF4B1PARAREQVQxc4AuHVxuEXuAyDl/AE56y/l/AEgW68fCWaFrH+uTg6S5mv1WXJVK5ejjdboS+TWVOkz94kqOXN/fmxQ3xYdEvueF9H3rQ3r16k2JK9Tta1/r3FLvg9ARLgPQEgorxsq441AHY/j63O1To4wFZoz3QdACeWAXbkccL+jz/P3C3n805Fqcbk9EX1p8J6w3Lc3pP9eDoSFe42E5feaCwmlhG15f1FSaV6xbK+hznAJdHuF65y9v6od2J6/a6P9g6tiBwBEREQmqsRrABSvASAiIvIfXgPgHa8BICIiCiAOAIiIyLccZZf88OKxxx5DQ0MD4vE4Fi5ciFdeeUV8fUdHBxYuXIh4PI5LLrkEjz/++Jj2N954A1/60pcwZ84cWJaFrVu3elqvyeAAgIiIfMuFBRd2CY/JnwJ47rnnsGbNGqxfvx6HDx/G9ddfj5tuuglvvvlm0dcfP34cN998M66//nocPnwYDz30EO6//348//zzhdcMDg7ikksuwSOPPIJUKuV5e0wGBwBERESTsGXLFnzjG9/AN7/5TVx++eXYunUrZs+eje3btxd9/eOPP46PfOQj2Lp1Ky6//HJ885vfxNe//nU8+uijhdd84hOfwA9+8APcfvvtiBmqVU6Vir0IMF8bBYrEbbL1JZT0FWJ+o+36tly9UNJXiPkBAOr0EbW4EPMDgESNPpI3rUofx5smxPwAYFp0QNt2vhDzA4Dzw/q+dbY+6ldjy2WGpZK/UswPkMv62kKsLlRCuV/T1KHrcXwtrS9QvnWW1td0oZQUE5RKBQPAQJG47xmZsD4iKEVOAaBGaK8Oy7HSKqH9vbA+ONZry5HgEVsfE8xbhsOxJe1P0udj+gtXiAmaSuVq2vP5s1djd6ouAuzr6xvzfCwWK/qLOJvN4tChQ/jud7875vnGxkbs3bu36M/Yt28fGhsbxzz32c9+Fk8++SRyuRwiEfk7Ui6cASAiIt+aqmsAZs+ejUQiUXi0tLQU/XmnTp2C4zhIJpNjnk8mk0in00X7pNPpoq/P5/M4derUFGwFbyp2BoCIiOhs6erqQn39X6aBTdPw1oduvqSUGvec6fXFnj+bOAAgIiLfGr0IsIRiQO/3ra+vHzMA0Jk+fTpCodC4v/a7u7vH/ZV/RiqVKvr6cDiMCy64wOOal46nAIiIyLfc928F7PUx2Wt1otEoFi5ciPb29jHPt7e3Y8mSJUX7LF68eNzr9+zZg2uuueacnf8HOAAgIiKalLVr1+LHP/4xfvKTn+Do0aP49re/jTfffBN33303AGDdunX4+7//+8Lr7777bpw4cQJr167F0aNH8ZOf/ARPPvkk/umf/qnwmmw2i87OTnR2diKbzeLPf/4zOjs78cc//rFs74OnAIiIyLdKuZnPaP/JJxZuu+02vPvuu/je976HkydPYv78+XjxxRdx8cUXAwBOnjw55p4ADQ0NePHFF/Htb38b27Ztw6xZs/CjH/0IX/rSlwqveeutt7BgwYLC/x999FE8+uijWLp0KV566SXP709SsQOAXE0YKjJ+9bK18gedqxXa6gw/U6jqJ0b9hJgfAFTV6aNx5wkxPwC4QIj6XRjr17ZNF9oAYHoko22bFtLH/ACgLqRfZynqZ6roJ7VLMT8AsIVonKmvRIoXuYaDjtdo0l/bex225SnOOqGSYJ2rvxBLipwCcvXJqpC8L8Zs/Xc6bAvbX2gDgNNC+xD0kUcAyEuHa+HzsQzFbqRUqeXKn7vlFF+nfO7s/WpxPUzjj+3v7TuzevVqrF69umjb008/Pe65pUuX4rXXXtMub86cOYULA8+Wih0AEBERmTjKEktZT6R/UPEaACIiogDiDAAREfnWmav5vfc/u9PulYQDACIi8i1X2cZrVOT+wR0A8BQAERFRAHEGgIiIfIunALzjAICIiHzLRWlX8nuvrel/FTsAyNXaUJHxo7pcraEsaZ2+PVcrj/S8Zv2lnD8AnC9k/WdU6/P4AJCq0rdfGNW3TQ/L9wGYJrSfZ7oPgJC/lrL8cUtfRhUwlfSVP7uQ5W0UbzpwSPcYdwxFPFyPByXb8F58t50gl7weVvpytOUqLx03lJeW+sZCwj0CDPt4KYEz6T4BeaEcszHL7woFbAy/HW2n+LLzOZ5d9oOKHQAQERGZlH4joOAOVjgAICIi3yr9VsDBHQAE950TEREFGGcAiIjIt1xY4jUoE+kfVBwAEBGRb/EUgHccABARkW+Vfh8ADgAqTq7aghsdPzWTqzHEAGv0bfkaOQKlavURnniNPspkKukrRf2kmB8AzIqd1rZNDwsRQaENkKN+pjKrNUKEKiLkhiLGiJrYXBamm4DkhOlB2xSNM8QEtcstIcpXtm3sMT4IAI6hqxQPzSn9vlZKeWlT34i4j8tRP69MW9gV4nzDQpQv7xjKAQvtVt7QV7OZnGxwp9X9pGIHAERERCausjzfc+NM/6DiAICIiHzLLfEUQJDvAxDcd05ERBRgnAEgIiLfKr0ccHD/DuYAgIiIfMuBBaeELH8pff0uuEMfIiKiAOMMABER+RZPAXhXsQOAfLUFFRs/NZMXcv6j/fRpWrdGzu9GqvVZ/4SQ9b+galBcrpT1l3L+ADAj0qdvC+vbpoXkcsD1lr6UarVtKJUq5tD1Qh5z8RPhKP06iZ+6YZVcYbnmvt7er6lkr1QO2JTzlw51+qK8JX52hq7SZ5cTMvcRJdeqjQqV3k1Z/pBwPwtTuWav8kJZZABwhPsAOEKWP2fI8ueFdtvQV9fuhM/etLqD0qbxy3NXB38I7tCHiIgowCp2BoCIiMiEpwC84wCAiIh8i8WAvOMAgIiIfEuVWA5YMQZIREREQcIZACIi8i2eAvBu0gOAl19+GT/4wQ9w6NAhnDx5Ert378att95aaFdKYcOGDdixYwd6enqwaNEibNu2DVdeeeWkfk6+GlCxYs/LERxHiPqFauR4W221Pho3TYj6JeNy6d1kVIjyCTE/AEiFT2vbLhCifufZ+vcCANVCDCpumBGLCJGwUJmm00xle71GeUylaqV4kWOI+XmNJpnKDEuLtaXYIgDb48djOkSW8rlL+5NU3thU+jhi60v+miKEISFCaAttJtIvm5wrxwDzQgwwm9f37XMM8UIh6pcXlgsAdk4TA5S7TSlWA/Ru0kOfgYEBXH311WhtbS3avnnzZmzZsgWtra04cOAAUqkUli9fjkxG/iVJREREZ8+kZwBuuukm3HTTTUXblFLYunUr1q9fjxUrVgAAnnnmGSSTSezcuROrVq0qbW2JiIg+wCmxHHApff1uSt/58ePHkU6n0djYWHguFoth6dKl2Lt3b9E+IyMj6OvrG/MgIiKaiDOnAEp5BNWUDgDS6TQAIJlMjnk+mUwW2j6spaUFiUSi8Jg9e/ZUrhIREREVUZa5D+tDF/UopcY9d8a6devQ29tbeHR1dZVjlYiI6K+QC7vkR1BNaQwwlUoBGJ0JmDlzZuH57u7ucbMCZ8RiMcRiRS73JyIiMnCUZUzkmPoH1ZQOABoaGpBKpdDe3o4FCxYAALLZLDo6OrBp06ZJLSsfB1R8/PNOlSG3VSXE26r01f4AIFE1rG2bFtPHAC+MygmH6RF9u1TRD/Ae9auz5WBctRS9suQRsdfIlynKJ1XeMwWvpHebE35szvBeclJsy5Af9hoDNN3VLCJtDaGKHQDAY3VDU6pLarcNlQSl/UmqLhkxVOWTon4hZQqOyjFaHdNfkzmh4t+IKx+Ohxz91hjK69uyeXm5A1mhyqDQBgD5rCYGGNzfqb4y6QFAf38//vjHPxb+f/z4cXR2dmLatGn4yEc+gjVr1mDjxo2YO3cu5s6di40bN6K6uhorV66c0hUnIiLifQC8m/QA4ODBg/ibv/mbwv/Xrl0LALjjjjvw9NNP48EHH8TQ0BBWr15duBHQnj17UFdXN3VrTUREBECVWA1Q8U6AE7ds2TIoYRrRsiw0Nzejubm5lPUiIiIycmB5Pt12pn9QBXfoQ0REFGAsBkRERL7lqtLO47uG68r/mnEAQEREvuWWeA1AKX39LrjvnIiIKMAqdgbAjSsgPn5uxhVy/gAQqdKXAK0Tcv4AcL6Q9Z8e0+fxp0f0bQBwoZD1nybk/AE5658Qsv41hiy/lPW3yzQudAxFe6X7BGQNZW6lrP+wlOU35ral+wDI6fisMT1fXNS4nfTtpguaHPGOClKb9zLDtmGdvO5vxhLFYtlq+X4JtngfDf130oH8fc6K9wGQ7noADLpRfZtwH4Bhw30Asln9OuUM9wFwRopvY9ewv0wlF5bx3hmm/kFVsQMAIiIiE94J0DueAiAiIgogzgAQEZFv8SJA7zgAICIi33JR4q2AA3wNQHCHPkRERAHGGQAiIvItVWIKQAV4BqBiBwBOTEEViQFacTkiVRXXxwDro3KJzwukGKAQ9ZsWNkX59Ms9z5ajiVJZXynqF7Pkj1aKSJUiJ5RZzQnlWUfb9dEhKeYHAINCvGpYaDNF+YaVUGbV0NdUGlb7Mw0RtagllLy29Ps/AESEvq4QL6w2RBNDYuxLfj/SvhixvEUpASAkxDDN0cO8tsUVvpMO5O/zsK3fnwbDMbGvFAMciOn7DuQMy43rl5sfkY8jTrz4NpZKe081VgP0rmIHAERERCa8CNC74L5zIiKiAOMMABER+RZPAXjHAQAREfkWbwXsHU8BEBERBRBnAIiIyLd4CsA7DgCIiMi3OADwrmIHACrmQsXG54cjMX0+FwCqY1ltWyI2JPY9L6LP658fHtC2GUv6hvTLrbEM70fISEtZ/1Ly0yZes/5Szh8AhoVmKecPAANKvy2kLP+woQSr1Dcn/MzRdm+fgZTVH23X7zNZw+cet/X3CXAg30NAJq2z/LlLpXmlvH4p+7hUvthE2k45w/dZOhZI+xoADAj3CeiNVOnbYnFxuX1Z/XKHYvI65WLFvwPKle/9QJWhYgcAREREJpwB8I4DACIi8i0OALxjCoCIiCiAOANARES+pVBalv/sVS2oPJwBICIi3zpzCqCUhxePPfYYGhoaEI/HsXDhQrzyyivi6zs6OrBw4ULE43FccsklePzxx8e95vnnn8cVV1yBWCyGK664Art37/a0bhPFAQAREfnWuRgAPPfcc1izZg3Wr1+Pw4cP4/rrr8dNN92EN998s+jrjx8/jptvvhnXX389Dh8+jIceegj3338/nn/++cJr9u3bh9tuuw1NTU347W9/i6amJnzlK1/Bq6++6nnbmFhKncW6jRPQ19eHRCKBi7Y/DLtqfHyl9jw5ypeqz2jb5tS+J/a9uOpdbVtD7B1t2/8J94jLTQoxwfNsOS5TZ+vP0sQsfUTHLmFKLG8o/Tqs9FGnESEGOODKu5oU9ZNifgCQcfVRJyleNeDKpVKlmKApBuh4HF9LsThAjgFKMT8AqLH1JbGlUsJ1hrLVUpy12hBrrBEyeTGh5HXcUPI6LJQDNnGFieERpd9OGVeOAZ529e/nbadW7Pvn/PnatuMjF2rbTgxdIC73T/3TtG3pvjqxb//p4vFDd2gY/+9bG9Db24v6+npxGV6d+V2x7P9+C+Ea+XssyQ+M4KXPb5/Uui5atAgf//jHsX379sJzl19+OW699Va0tLSMe/13vvMdvPDCCzh69Gjhubvvvhu//e1vsW/fPgDAbbfdhr6+Pvzyl78svObGG2/E+eefj5/97Gde356IMwBERORbZ3sGIJvN4tChQ2hsbBzzfGNjI/bu3Vu0z759+8a9/rOf/SwOHjyIXC4nvka3zKnAiwCJiMi3pioG2NfXN+b5WCyGWGz8zMKpU6fgOA6SyeSY55PJJNLpdNGfkU6ni74+n8/j1KlTmDlzpvY1umVOBc4AEBFR4M2ePRuJRKLwKDaV/0HWh+7SqpQa95zp9R9+frLLLBVnAIiIyLeUsqBKmAE407erq2vMNQDF/voHgOnTpyMUCo37y7y7u3vcX/BnpFKpoq8Ph8O44IILxNfoljkVOANARES+5cIq+QEA9fX1Yx66AUA0GsXChQvR3t4+5vn29nYsWbKkaJ/FixePe/2ePXtwzTXXIBKJiK/RLXMqcAaAiIhoEtauXYumpiZcc801WLx4MXbs2IE333wTd999NwBg3bp1+POf/4z/+I//ADB6xX9rayvWrl2Lu+66C/v27cOTTz455ur+Bx54ADfccAM2bdqEW265Bb/4xS/w3//93/j1r39dtvdRsQOAUNSBHRsfH6qKyjGnuog+rlQfliOEiZC+vc7Wt1UL0SoAiFv6WFdciDkBcsWzUqJ+UszJMSRDpYp/w0LfUir6STG/0XZ9NTQp6meKAUoV/UYMlQQdj59PyHBvspgQ9RtWUbGvFF2UIoJGhpigJCRUl7SFSGTEsJ1sS99u+u5I7dJ3Uvqum9pNxxHpGCQdu0zHPemY2RuVv3dDseL7m+XI0c+pdC5qAdx2221499138b3vfQ8nT57E/Pnz8eKLL+Liiy8GAJw8eXLMPQEaGhrw4osv4tvf/ja2bduGWbNm4Uc/+hG+9KUvFV6zZMkS7Nq1C//8z/+Mf/mXf8Gll16K5557DosWLfL83kwqdgBARERkMlXXAEzW6tWrsXr16qJtTz/99Ljnli5ditdee01c5pe//GV8+ctf9rQ+XvAaACIiogDiDAAREfkWywF7xwEAERH51rk6BfDXgAMAIiLyLVXiDECQBwC8BoCIiCiAOANARES+pQCUUtO2osrhnmUVOwAIhR2EwuOzpPGwXG6zOqzPSNeGDTlbIUsrZaSlUqgAEBdmmCKG+wBIJU1DQl9HyOoDgCvkq3OGcsA54ds2LEynDRvuAyBl/aWcPwD0Od76mrL8UjngEVM5YOWxHLAhSx4T1slUDli6r4FTpqnQkJBfB4CIsK9KWf+cocxwSMjy24ZSwdJ3S/qNETGsk3QfANNxRDoGSccu03FPOmaajrfhIsdoAHA0z5eDCwtWSfdE4SkAIiIiCpCKnQEgIiIyYQrAOw4AiIjIt1xlweJ9ADzhKQAiIqIA4gwAERH5llIlpgACHAPgAICIiHyL1wB4V7EDgEjERShSpBywEFkBgBoh8lJtZ8W+UjnOuCVEZYzRH/2ZFqm0KGCII5VAKvkrlfsFgGFhxCxF/QYMpWql0rxSzA+Qo379Qt9BV14nMQboyl8fr+cWpTK2ABCz9dEsUzSx2hDF9Eos22v4fsgxQP17lfqN9tVvx3AJx3zpOxkxxAvl0uCm44j+GCQdu0zHPemYaTreRoocowHAjsifDVWGih0AEBERmXAGwDsOAIiIyLeYAvCOAwAiIvItXgToHWOAREREAcQZACIi8q3RGYBSrgGYwpXxGQ4AiIjIt3gRoHc8BUBERBRAFTsDEA3nEQqPz9TGDOUpq0JSVlbOw8pZf/3PjRkGkFLJX7tMYzDXUOVaLAdsmBMbFsrcDgq5+UEh529qN5UDlrL+GaFtyJHLAQ8J9wnIuXLmu1z3Acja+q9t3pb3JydUnv0tKmTYpTYAiArlpyO2kJs3lbwWMvem74fXuyWYvs/SsSBmKAMtHYPkewTIxz3pmGk63kY17Y6h31RSECs0T6h/UFXsAICIiMiEpwC84ykAIiKiAOIMABER+RfPAXjGAQAREflXiacAEOBTABwAEBGRb/FOgN7xGgAiIqIAmtQMQEtLC37+85/jD3/4A6qqqrBkyRJs2rQJ8+bNK7xGKYUNGzZgx44d6OnpwaJFi7Bt2zZceeWVk1qxiO0iFBofi4kLkRVALpUas+W+Xkv+miJDUjTIRnmmn6SYHyCX/JW3EpAT3s+w0sfqpHK/pnZT2V4p6tfv6JdrjAEK7XlTDNDjZ2sbTkqGbf2+6ITkn+l4HPPbxoiafhubygGLETZhb8yVsI9HDO/HaxDQ9H2WjgWmnygdg6RtaDruScdM0/E2UuQYDQC2EN+cakwBeDepo0FHRwfuuece7N+/H+3t7cjn82hsbMTAwEDhNZs3b8aWLVvQ2tqKAwcOIJVKYfny5chkMlO+8kREFHDKKv0RUJOaAWhraxvz/6eeegozZszAoUOHcMMNN0Apha1bt2L9+vVYsWIFAOCZZ55BMpnEzp07sWrVqqlbcyIiIvKspGsAent7AQDTpk0DABw/fhzpdBqNjY2F18RiMSxduhR79+4tuoyRkRH09fWNeRAREU3EmYsAS3kElecBgFIKa9euxSc/+UnMnz8fAJBOpwEAyWRyzGuTyWSh7cNaWlqQSCQKj9mzZ3tdJSIiCho1BY+A8jwAuPfee/H666/jZz/72bg2yxp7TkUpNe65M9atW4fe3t7Co6ury+sqERER0QR5ug/AfffdhxdeeAEvv/wyLrroosLzqVQKwOhMwMyZMwvPd3d3j5sVOCMWiyEWk68OJyIiKoYpAO8mNQBQSuG+++7D7t278dJLL6GhoWFMe0NDA1KpFNrb27FgwQIAQDabRUdHBzZt2jSpFQvZLsJFoiRhQ3xHihxFhWpao+36vlI0K6KZ3TgjJLSHhOpgJo4Qc3IMJ7Yc4f3kDFNiOaEaYFbpw0w5oQ2QqwGOCFUGR9v1u7IU5TPFAIfFGKD82bnCdpKYIndhj8s1kb47pu0/KMQAa+wRsa+8z+jfq2k/dYSqiubvh/4zkL6zpu9zSIgJmo4j0jFIrsYoH/ekz910vC12jAYAnMUYIIBAT+OXYlIDgHvuuQc7d+7EL37xC9TV1RXO6ycSCVRVVcGyLKxZswYbN27E3LlzMXfuXGzcuBHV1dVYuXJlWd4AERERTd6kBgDbt28HACxbtmzM80899RTuvPNOAMCDDz6IoaEhrF69unAjoD179qCurm5KVpiIiOgMngLwbtKnAEwsy0JzczOam5u9rhMREdHEsBqgZywGREREPma9/yilfzCxGBAREVEAcQaAiIj8i6cAPOMAgIiI/IsDAM8qdgAQ1t0HQCiFCsiZVlNZUil/HRXaQoYzKVIJ0HNFSuk6hnNiUp4/JxQ1lUoFm5Y7XMJ9ALIe2wAg6+jb84Y8vuvx6mLbkCV3hfs/2EL2HQBCQru0DU3bv9rO6vuaPndhn5H2CdN+epaT6BMilwOW3490DJKOXabjnngfAMPxtmLuA0CeVOwAgIiIyKjUkr6MARIREflPqRX9WA2QiIiIAoUzAERE5F+8CNAzDgCIiMi/eA2AZzwFQEREFEAVOwNgWQpWkchSpIRywCFDMMjU7ieu4b1I5VAdQ7xNil/llH6XMi1XjBcaSgnnpb6uvi3rmJarX2dTOWDT+9UJGfZxadhuej9SeVdxGxq2vxjXM37uwj4j7GuO4S83qVy2a9rGPiMdu0zHNTk6LfctdoyWni8HS40+SukfVBU7ACAiIjLiNQCecQBARET+xWsAPOM1AERERAHEGQAiIvIvngLwjAMAIiLyLw4APOMpACIiogDiDAAREfkXZwA8q9gBgG0pY2nTqRYS9gRpqsS2zs1VpG6Z9lzXUJbUEbaGlPmW+pn6mkrrSplwqa9bQklfZVgnz+WATWVuy/R+vG7D0b7l+dylvqb9tBTSd0u+I4J3puOItBWlY1c56Y7RZ/XYzRSAZzwFQEREFEAVOwNARERkwjsBescBABER+RevAfCMpwCIiIgCiAMAIiKiAOIpACIi8i0LJV4DMGVr4j+cASAiIv86EwMs5VFGPT09aGpqQiKRQCKRQFNTE06fPi2/JaXQ3NyMWbNmoaqqCsuWLcMbb7wx5jU7duzAsmXLUF9fD8uyjMsshgMAIiKiMlm5ciU6OzvR1taGtrY2dHZ2oqmpSeyzefNmbNmyBa2trThw4ABSqRSWL1+OTCZTeM3g4CBuvPFGPPTQQ57XjacAiIjIvyo4BXD06FG0tbVh//79WLRoEQDgiSeewOLFi3Hs2DHMmzdv/Oooha1bt2L9+vVYsWIFAOCZZ55BMpnEzp07sWrVKgDAmjVrAAAvvfSS5/XjDAAREfmXmoJHmezbtw+JRKLwyx8ArrvuOiQSCezdu7don+PHjyOdTqOxsbHwXCwWw9KlS7V9vOIMABERBV5fX9+Y/8diMcRisZKWmU6nMWPGjHHPz5gxA+l0WtsHAJLJ5Jjnk8kkTpw4UdL6fBhnAIiIyLfO3AmwlAcAzJ49u3ChXiKRQEtLi/ZnNjc3w7Is8XHw4MHR9StS40EpVfT5Me/rQ+0T6TNZnAEgIiL/mqJrALq6ulBfX194Wvrr/95778Xtt98uLnbOnDl4/fXX8fbbb49re+edd8b9hX9GKpUCMDoTMHPmzMLz3d3d2j5eVewAwFWW52pqXjlCItQV+rnKsPeV6W2YqsZ5X678fkLC1ghZQpu4FeW+pupiIaFd6msLP3O0XahGZ1gn03bUsUzLLen9CFXuPG7D0b7l+dylvl6370SU67slMR1HpK0oHbvKSXeMPtvH7qlQX18/ZgAgmT59OqZPn2583eLFi9Hb24vf/OY3uPbaawEAr776Knp7e7FkyZKifRoaGpBKpdDe3o4FCxYAALLZLDo6OrBp06YJvpuJ4SkAIiLyrwq+CPDyyy/HjTfeiLvuugv79+/H/v37cdddd+Hzn//8mATAZZddht27dwMYnfpfs2YNNm7ciN27d+N3v/sd7rzzTlRXV2PlypWFPul0Gp2dnfjjH/8IADhy5Ag6Ozvx3nvvTXj9KnYGgIiIyKTSqwH+9Kc/xf3331+4qv+LX/wiWltbx7zm2LFj6O3tLfz/wQcfxNDQEFavXo2enh4sWrQIe/bsQV1dXeE1jz/+ODZs2FD4/w033AAAeOqpp3DnnXdOaN04ACAiIiqTadOm4dlnnxVfoz50+seyLDQ3N6O5uVnbx9Q+ERwAEBGRf5V6O18fXq8wVTgAICIi/6rgOwFWOg4AiIjItyr9GoBKxhQAERFRAFXsDIBSFlSRczM5JY9ZciqkbXMM4x1Tu5/YhvcSEu4oJeXBASAkzJlFrLywXDkPHrEcT20AEJb62vq2aEjeTq60vxl2F9t0fwhdP8P2DwvbMRqSt5O0LcRtaNj+Urv5cxf2GWFfM+6nwj5u+n74jXTsMh3XpGOm6Xhb7BgtPV8WPAXgWcUOAIiIiIxKPAUQ5AHAX9cwmIiIiCaEMwBERORfPAXgGQcARETkXxwAeMZTAERERAHEGQAiIvIt3gfAu4odAORdG3DHT1DkXX1kBTBFWuS+UuQrK7Q5hjkkVyzkKa9TuUhTP1L0CjDE9aBvi1s5z8uN23LfEVu/K4/Y+piZY4gruSEhSlZk/xzT1xCh0i7XEJsL20IMUHivpvaY0Gba/uJnZ/rchX1GjBca9tNzUdLXRDoWmI4j0jFI2tdMxz2p3XS8zWu+A7rnqbLwUyIiIgqgip0BICIiMuJFgJ5xAEBERL7FawC84wCAiIj8LcC/xEvBawCIiIgCiDMARETkX7wGwDMOAIiIyLd4DYB3FTsAcHT3ASihHHBWyW83K/R1hUxxTsm5bUcoC+uI9wgAQpb+/Ypthgx0SMi/RwzfiIhUjraEkrLV9oi2bdjw2cXsiLatKiTn0L3KW4b7SnjModuGP0nCQklf03uV2qX7AMQM9wGQPjvT5y7vM/p9LWLYvNJ3QCoVPNru7eyoU8KxIGcoH+0KZ2ylY5fpuCfeB8BwvNXl/R3eB8AXKnYAQEREZMRTAJ5xAEBERL7FUwDecZ6GiIgogDgDQERE/sVTAJ5xAEBERP7FAYBnPAVAREQUQJOaAdi+fTu2b9+OP/3pTwCAK6+8Ev/6r/+Km266CQCglMKGDRuwY8cO9PT0YNGiRdi2bRuuvPLKSa9YzrXhOuPHJ8OOPu4FACOuUBbWlfsOK337sBCVcYRypoBcAtQ1DD+9Fgu2DWO7iBBzihiiTBHh/UilX2uEqBggx5GkzwYA3JC3sWzIGMMUYluGUqmuodSwjm2KYUoxQDsr9pVigHWhYW1btWG50mdr+tylfUba1+Q9Qt7HTd8Pr0zfZ7kcsEw6BknfD9NxTzpmmo63uSLHaODsxgB5EaB3k/qULrroIjzyyCM4ePAgDh48iE996lO45ZZb8MYbbwAANm/ejC1btqC1tRUHDhxAKpXC8uXLkclkyrLyREQUcGoKHgE1qQHAF77wBdx888342Mc+ho997GP4/ve/j9raWuzfvx9KKWzduhXr16/HihUrMH/+fDzzzDMYHBzEzp07y7X+REQUZBwAeOZ5nsZxHOzatQsDAwNYvHgxjh8/jnQ6jcbGxsJrYrEYli5dir1792qXMzIygr6+vjEPIiIiKq9JDwCOHDmC2tpaxGIx3H333di9ezeuuOIKpNNpAEAymRzz+mQyWWgrpqWlBYlEovCYPXv2ZFeJiIgC6sw1AKU8gmrSA4B58+ahs7MT+/fvx7e+9S3ccccd+P3vf19otz50j22l1LjnPmjdunXo7e0tPLq6uia7SkREFFQ8BeDZpO8DEI1G8dGPfhQAcM011+DAgQP44Q9/iO985zsAgHQ6jZkzZxZe393dPW5W4INisRhisdhkV4OIiIhKUPKNgJRSGBkZQUNDA1KpFNrb27FgwQIAQDabRUdHBzZt2jTp5WbzYYTy41dvpMhzHzQkxFYG3ajYV44BCvFCJUekpGqBUrWzUd6CgLahEp0Ug4oYKqXFhXXOCVXjspDjYHJFM6+BSFnIsP3Drr5dik8B5YsBSlX74saqffp9tVaIAdbZQ4bl6j9bqW20Xb/O0r5m2k+lfdz0/fBKivkB8rFgxPCXqHQMko5dpuOedMw0HW+zmnYnbwo1Th3GAL2b1ADgoYcewk033YTZs2cjk8lg165deOmll9DW1gbLsrBmzRps3LgRc+fOxdy5c7Fx40ZUV1dj5cqV5Vp/IiIKMt4J0LNJDQDefvttNDU14eTJk0gkErjqqqvQ1taG5cuXAwAefPBBDA0NYfXq1YUbAe3Zswd1dXVlWXkiIiLyZlIDgCeffFJstywLzc3NaG5uLmWdiIiIJoYzAJ6xGBAREfmW9f6jlP5BxWJAREREAcQZACIi8i+eAvCMAwAiIvItxgC9q9gBQC5nw82Nz34P5eXylAN5/U2FTHnYQVff12up4NF2fSY2IrQBQFi4D0BIKHdqEhIy1BHDmaG4pV/nnFDUtMYy3C/B1r9Xt0xnqyLCewGAmKXP3I/Y8tfHUR5LFBvuTSCtk+k+ADGhXcr61wv3CAAM5YANn7u0P8WFE7RSuV9A3sdL4QhZ/pzh+zws9DUfR6Ssv/fjnnTMNB1vc0WO0QDg5M7i2WXOAHjGawCIiIgCqGJnAIiIiCYkwH/Fl4IDACIi8i1eA+AdTwEQEREFEGcAiIjIv3gRoGccABARkW/xFIB3FTsAcPIhqPz4iMmwoTzloBBb6RfiLgCQCVdp2wZC+r4DthyRqhaiQXEhFgQYYmol7LhiOWBDCWKphHFcWKmcIXJXZ9iOXoWEb/iAJe8TUmxuxJUjUo7Hm4yGDB+stE6mWKMY1xPaTOWApc9OivmNtuvfr1Ty17SfSvu4iRT1ywtRV6ncLwAMCx/tgFDuFwAGhKhfxtEfu0zHPemYaTre5oscowHA1TxPlaViBwBERERGPAXgGQcARETkWzwF4B1TAERERAHEGQAiIvIvngLwjAMAIiLyLw4APOMAgIiIfIvXAHjHawCIiIgCqGJnAJxsCCpUpBxwVs5eZ3JxbVtfXp+VBYBeIUubcfVt9a6cXx+29LltqTwoIJcLtoVyqLYhgy61m8qoSuWC49C/H8eQBxeVcI+AiFA+V2oDgGEh6x8XPlcAcDyOr0PCNgTkdTaVA5ay/tL7Md2joUZYp2rjfQCErL+wj5v2U9N3QOJK97MQvpOm7/OwUCJaKukLyMcg6dhlOu5Jx0zT8dYZ0dwHIHsW7wPAUwCeVewAgIiIyMRSCpby/lu8lL5+x1MAREREAcQZACIi8i+eAvCMAwAiIvItpgC84ykAIiKiAOIMABER+RdPAXhWuQOAkRBgj4+SjIzIq9yfjWrbeoW4CwD0RvRxmdPham2bqVRqtRC9iip9G2CIAUolfS3vMRxTGVWpDKsUn5LKvo4qT0xQ2oYRw88ctvQxqJyhfGtOefsMTCV9xRigIZooxQSlvlLMD5CjfnFDGi8ixPWkfc20n0r7ookc9dO3DRp+ZEbp96c+Vz4+nXb0x6BeIepnOu5Jx0zT8RaaGKD2+TLgKQDveAqAiIgogDgAICIi/1JT8Cijnp4eNDU1IZFIIJFIoKmpCadPnxb7KKXQ3NyMWbNmoaqqCsuWLcMbb7xRaH/vvfdw3333Yd68eaiursZHPvIR3H///ejt7Z3UunEAQEREvnXmFEApj3JauXIlOjs70dbWhra2NnR2dqKpqUnss3nzZmzZsgWtra04cOAAUqkUli9fjkwmAwB466238NZbb+HRRx/FkSNH8PTTT6OtrQ3f+MY3JrVulXsNABERkUkFXwR49OhRtLW1Yf/+/Vi0aBEA4IknnsDixYtx7NgxzJs3b/zqKIWtW7di/fr1WLFiBQDgmWeeQTKZxM6dO7Fq1SrMnz8fzz//fKHPpZdeiu9///v4u7/7O+TzeYTDE/vVzhkAIiKiMti3bx8SiUThlz8AXHfddUgkEti7d2/RPsePH0c6nUZjY2PhuVgshqVLl2r7AEBvby/q6+sn/Msf4AwAERH53FRM4/f19Y35fywWQywmF2gySafTmDFjxrjnZ8yYgXQ6re0DAMlkcszzyWQSJ06cKNrn3Xffxb/9279h1apVk1q/ih0AWCM2LHv8BEXeEEsZHBFigFG5KtbpiD5m0xOu0bZJFdYAOV4VN0W+lD5+FYIczZKYKql57StVCnSMc21Su/eIYEhYbtSWlxsVIl9ZQ8zP9TjBZhuqAUbFyJ0cA5QihtK+aK7oJ/1MQ3VJjxX/XMN2kjiGAjAjwvduUOg7YIiGSlG+95xasW9PXn8MOp0TIoIj8nFPOmaajrfWSPHPTvd8WSg1+iilP4DZs2ePefrhhx9Gc3Nz0S7Nzc3YsGGDuNgDBw4AAKwi+7BSqujzH/Thdl2fvr4+fO5zn8MVV1yBhx9+WFzmh1XsAICIiOhs6erqQn19feH/0l//9957L26//XZxeXPmzMHrr7+Ot99+e1zbO++8M+4v/DNSqRSA0ZmAmTNnFp7v7u4e1yeTyeDGG29EbW0tdu/ejUhELt/8YRwAEBGRb03VjYDq6+vHDAAk06dPx/Tp042vW7x4MXp7e/Gb3/wG1157LQDg1VdfRW9vL5YsWVK0T0NDA1KpFNrb27FgwQIAQDabRUdHBzZt2lR4XV9fHz772c8iFovhhRdeQDwu3/CpGF4ESERE/lXB9wG4/PLLceONN+Kuu+7C/v37sX//ftx11134/Oc/PyYBcNlll2H37t0ARqf+16xZg40bN2L37t343e9+hzvvvBPV1dVYuXIlgNG//BsbGzEwMIAnn3wSfX19SKfTSKfTcJyJny7lDAAREVGZ/PSnP8X9999fuKr/i1/8IlpbW8e85tixY2Nu4vPggw9iaGgIq1evRk9PDxYtWoQ9e/agrq4OAHDo0CG8+uqrAICPfvSjY5Z1/PhxzJkzZ0LrxgEAERH5luWOPkrpX07Tpk3Ds88+K75GfegiRsuy0NzcrL0IcdmyZeP6eMEBABER+VcF3wio0vEaACIiogCq2BmA0IgFu0jmMT8sZ6+HYkK5TcNNHd4d0Wdpa8L6rH+1nRWXK5UDjrtybjsEfalhW8ywy/cIkPL6pjKrErlEsWmorZ+LM91DoFq4T4At9M0Z5v+kUsKmcr9ZoZStJGq454GU5Y+Y3o+wjeNC34ippK90bwgh5w/I+4wpry+R7hOQU/J2GhDaM67+cz1tKunrCvcByMv3ATiV07dLx66+rHzcGxrWHzOV4XgbHi7+uVsj3u8zMlksB+xdxQ4AiIiIjKboRkBBxAEAERH5FmcAvOM1AERERAHEGQAiIvIvpgA84wCAiIh8i6cAvOMpACIiogCq2BkAe9hCCOOjJO6QHEvJRfWRlkxUjujEw/roXHVYH9czxQBjtr6vVNoVAEJiyVN9vNA1lLmtFua9THG9Yp/LVJCWGzX+SGmd9dtQKhVsajf1jXgsYWxcrhjXM8UApc9d3y9qKF1ayj7htayvKRoqRf2kkr6AKeqnj9WZSvq+k9cXm5FifgBwakTf3iPEADND8nEvN6Q/ZtqG421IEwOE7vlyYArAs4odABAREZnwFIB3PAVAREQUQJwBICIi/2IKwDMOAIiIyLd4CsA7ngIgIiIKIM4AEBGRf7lq9FFK/4DiAICIiPyL1wB4VrEDgPAwECrywbiGQLgb1edWh6NRsW9vWJ+XrRLuA1AVkkv6xmz9/QWk0q4AYHvMSEv3CAAAR/i5ccM6SaVfy3WPANO5KjGtLK2SKQMs5OqlMsMAkPN4hk0q2QsAIeGkpZTzB+Ssv7QNSzlX6Bq2sSuss5T1zxmWOyw0DxpKOUtZ/3eFrH+3kPMHgFO5Om3bO1l9GwC8J2T9e4Ws//CQfNyDkPUPDcnf57CmWrklH36mlIUSrwGYsjXxH14DQEREFEAVOwNARERkxDsBesYBABER+RZjgN7xFAAREVEAcQaAiIj8iykAzzgAICIi37KUglXCefxS+vpdSQOAlpYWPPTQQ3jggQewdetWAIBSChs2bMCOHTvQ09ODRYsWYdu2bbjyyisnt2KDQKhIGs2VckwA3LA+0uKE5Q+6PySU+Qzro3Fh21CCVSjNGzKUb5W4whkcB/1i35yQ08kJsUUAiEslZaXyuYaSsudCyLBKYoStlL4CKeYHGEoUl2kTm/ZSKcpn4gjbSQrYDiv5zQ66+sNbn9J/1wG5rK8U9evOyTHAt7P69reHDTHAIX0MsH9Q/36cAfkwHxrQHzPDg4YY4GDx589mDJC883wNwIEDB7Bjxw5cddVVY57fvHkztmzZgtbWVhw4cACpVArLly9HJpMpeWWJiIjGcKfgEVCeBgD9/f342te+hieeeALnn39+4XmlFLZu3Yr169djxYoVmD9/Pp555hkMDg5i586dU7bSREREwF9OAZTyCCpPA4B77rkHn/vc5/CZz3xmzPPHjx9HOp1GY2Nj4blYLIalS5di7969RZc1MjKCvr6+MQ8iIiIqr0lfA7Br1y689tprOHDgwLi2dDoNAEgmk2OeTyaTOHHiRNHltbS0YMOGDZNdDSIiIqYASjCpGYCuri488MADePbZZxGP6+89bX3oYi+l1Ljnzli3bh16e3sLj66ursmsEhERBdmZOwGW8gioSc0AHDp0CN3d3Vi4cGHhOcdx8PLLL6O1tRXHjh0DMDoTMHPmzMJruru7x80KnBGLxRCLyVfkEhERFcM7AXo3qQHApz/9aRw5cmTMc//wD/+Ayy67DN/5zndwySWXIJVKob29HQsWLAAAZLNZdHR0YNOmTZNbsUGFUH78J+OGTTFAfbsSIoIAkAvrq2b1FitN+L6QIQYYNlTX88oVYlCOkid3cqEBbVsWw2LfGkuqbqjfFhHDSLtcETaJY/jyO0LWzzHE0LwyLldoto3VDfVN0l5sqrwnMW3jnLBSOWE/HlDy4Svj6mcpTzs1Yt938vpI3imh7a2R88Tlpof0fd8VYn4A0DtQpW3LDeqPXbYQ8wPkqF9Yf5h4v734h2tlA/xb1UcmNQCoq6vD/PnzxzxXU1ODCy64oPD8mjVrsHHjRsydOxdz587Fxo0bUV1djZUrV07dWhMREQEsBlSCKb8T4IMPPoihoSGsXr26cCOgPXv2oK5OvskFERHRZFnu6KOU/kFV8gDgpZdeGvN/y7LQ3NyM5ubmUhdNREREZcJaAERE5F88BeAZBwBERORfvA+AZ55rARAREZF/cQaAiIh8i+WAvavYAUBkUCGcG//BGKK/UNJ9AAxB81xIn5cdsfU529OG+wCUK96eU/r1zRk21LCK6NuEewQAwLCtv09A3NIXcI0b7ocQEr6ItmGezlRCV8eUuXel+wAYPlnpPg0S2/BepDLDjmk7eZzvLOd2Ghb2Y2k/lXL+gJz1fy+vL/cLAKeE9ney+kSTlPMHgO5BfftpIecPACMD+mOQ1S+U9B2Qt39E+LpHNDn/M6Ka9nyRY3fZ8BoAz3gKgIiIKIAqdgaAiIjISEG+jeVE+gcUBwBERORbvAbAOw4AiIjIvxRKvAZgytbEd3gNABERUQBxBoCIiPyLKQDPKnYAEOl3EY6Mv7LDDcmTFlIVXGXLcRhl6zvnLf2mGoIcRxJ/pqE9L0SkRhz9Og3H9PEpABgWYoLDrtx3IDSkbauxR7RtUkTQ1G6Kr9lCRQ+v0TdAjrC5hpLLpvibjhSHBEp8r0Jzud6rFOUztQ+4MW1bxpFjc+8JMcBTOTmud2pEiAEKbaaSvlLUbyhjOI5k9N/ZcL/+84n0G2KAGf1OEemX96dIf/F90cqdxQo7LkrLWge4GBBPARAREQVQxc4AEBERmTAF4B0HAERE5F+8BsAzngIgIiIKIM4AEBGRf3EGwDMOAIiIyL84APCMpwCIiIgCqGJnACIDeYTD+XHPq5C8ysrW5+YNUWYoSwiTWsI9AgybUbpPgOsa8tVCe15oG3EN6+Tos9f9ETmPfH5YXz+0TigVLN0jAADitv4+ABFr/L7wQVL+3RaCviEhU2/iGHYo1+P4WlpfoHzrLK2v6Z4GUvlp430lpKy/UPK3J6/P+QNAT06fyX8vK/d9b0ToK2T9e0so6Svl/AFD1j+j/3wiGXGxiPTr26KanP9f+hb/Xlp5+fs6pXgfAM8qdgBARERkwhigdxwAEBGRf/EaAM94DQAREVEAcQaAiIj8y1WAVcJf8W5wZwA4ACAiIv/iKQDPeAqAiIgogCp2BiDcn0U4VCTbYYp7iO36iKC5s9BmiIPlxYiU/IZcoT2b17+fobwcvRoU2gei+lgWAPRH9O21YX3UrzakjwgCQLWd1baZSglLMcGQMD0YKmMGyPE4vj5X6+Qo/b4mxfwAuaTvoCtE3wD0O/qoX39ev6+dFmJ+ANCT1UfyeoSYHwD0DunXaWBIv05ZKeYHwOrXf2elmB8gR/2iffp+0T75L9xYxhH6ynG+cEbzfXf03+WpV+IMQAnlwv2uYgcARERERjwF4BlPARAREQUQZwCIiMi/XIWSpvEDnALgDAAREfmXckt/lFFPTw+ampqQSCSQSCTQ1NSE06dPy29JKTQ3N2PWrFmoqqrCsmXL8MYbb4x5zapVq3DppZeiqqoKF154IW655Rb84Q9/mNS6cQBARERUJitXrkRnZyfa2trQ1taGzs5ONDU1iX02b96MLVu2oLW1FQcOHEAqlcLy5cuRyfylsMPChQvx1FNP4ejRo/iv//ovKKXQ2NgIx9Ff1PlhPAVARET+VcEXAR49ehRtbW3Yv38/Fi1aBAB44oknsHjxYhw7dgzz5s0rsjoKW7duxfr167FixQoAwDPPPINkMomdO3di1apVAIB//Md/LPSZM2cO/v3f/x1XX301/vSnP+HSSy+d0PpV7ADAzgzDDo3/YMq7wkJMUIhIWUIbAFhSRT9H7pvN69sdIQaYc+TI43BevyUHcnIMsDeqj0jVRfQxwLqwHAOsCumjflKlQACICTHAiC1EBA3nDu0SKu+dC64hkipV9csJFSRHTDFAoeKfVHkSADJ5/f6UEfbFTFauWtmX1fftH5b38aFBfbszoN8W9oD8vQsPCFX7+uVjQUSK+mX0+7EU8wOAaK/++xHukyt42pni32nbkftNqQq+BmDfvn1IJBKFX/4AcN111yGRSGDv3r1FBwDHjx9HOp1GY2Nj4blYLIalS5di7969hQHABw0MDOCpp55CQ0MDZs+ePeH14ykAIiLyrzMzAKU8APT19Y15jIyUPohJp9OYMWPGuOdnzJiBdDqt7QMAyWRyzPPJZHJcn8ceewy1tbWora1FW1sb2tvbEY3K96L4IA4AiIgo8GbPnl24UC+RSKClpUX72ubmZliWJT4OHjwIALCs8TM7Sqmiz3/Qh9uL9fna176Gw4cPo6OjA3PnzsVXvvIVDA/LM60fVLGnAIiIiIwUSrwGYPSfrq4u1NfXF56OxfSnge69917cfvvt4mLnzJmD119/HW+//fa4tnfeeWfcX/hnpFIpAKMzATNnziw8393dPa7PmcHK3Llzcd111+H888/H7t278dWvflVctzM4ACAiIv+aoosA6+vrxwwAJNOnT8f06dONr1u8eDF6e3vxm9/8Btdeey0A4NVXX0Vvby+WLFlStE9DQwNSqRTa29uxYMECAEA2m0VHRwc2bdpkeCtqUqcueAqAiIioDC6//HLceOONuOuuu7B//37s378fd911Fz7/+c+PuQDwsssuw+7duwGMTv2vWbMGGzduxO7du/G73/0Od955J6qrq7Fy5UoAwP/+7/+ipaUFhw4dwptvvol9+/bhK1/5CqqqqnDzzTdPeP04A0BERP7lukApBbTc8iZ9fvrTn+L+++8vXNX/xS9+Ea2trWNec+zYMfT29hb+/+CDD2JoaAirV69GT08PFi1ahD179qCurg4AEI/H8corr2Dr1q3o6elBMpnEDTfcgL179xa96FCHAwAiIvKvCr4PAABMmzYNzz77rGEVxq6DZVlobm5Gc3Nz0dfPmjULL774YsnrVrEDACszAKtI9ts2fFjiGzJ9zkK7lOW3DCV9LSHrbws5fwDIifcB0LcN5OSzOyMj+i01GDeUbxXaeyP6vH5NWC4RWi20S/cIAIBYSMgyW/ocdERoAwBbKCUstQGAa7g/RDmWa/qZOaXPqeeFthFHPlRIWf/BvLw/DQjtgzn9co1Z/mF937zQBgBqUL8tQkLWPyLk/AEgPKBviwhZfkDO+kf79X/FSjl/QCjpC8DuGxL7Wpnib8hyz2Y5YPKqYgcARERERhU+A1DJOAAgIiL/quA7AVY6pgCIiIgCiDMARETkW0q5UCWU9C2lr99xAEBERP6lVGnT+LwGgIiIyIdUidcAcABQedz+frjW+HiQZfiwpJhgxJH7Wo4+VmQL5XMtVy4BagtJM1MM0BLa80I5YCcrX96RqxJKFGcNpV9H9BGqWEwo6RuR40jxsL69KmyIAQp9pRhg2Jan/8JCOWBTXK9cpKhf3lAOOC+VppZigML+DwBDef0+IZWeBoDhnFCGWNjXckKUFQDUkP792EIbAISH9NtYLOkrxPwAINIvHJ+ENgCI9Qlx1j5h/xdifgBg9QtRP03M7wy3v3i7qxgD9IOKHQAQEREZuS4gDNSNeA0AERGRD/EUgGeMARIREQUQZwCIiMi3lOtClXAKgDFAIiIiP+IpAM94CoCIiCiAOANARET+5SqglEhugGcAKnYA4A4MwrWKlAN25PKtltBuO4bMt6tvt9y48DPFxcLO6ydabEPZXjsnlBIW2pwR+f4CeaHdGZHXKRvTZ6hzcf0uNRSRt38kqs8yR8LyRo4K7eGQkJ823AdAyvpbhoOO1/sEmEr6qlLKAUv3AXD0n2tWuOcEAOSE9pzhvhKO8B1Qw/rlWob9NDykb5dy/gAQEqLxUtbflOWPDAhlezOG70dGfy+MUGZY22Zn5JK+KtOvb9Pk/M9wh4v/XFfJ9+2YUkoBKCUGGNwBAE8BEBERBVDFzgAQERGZKFdBlXAKQAV4BoADACIi8i/lorRTAIwBEhER+Q5nALzjNQBEREQBVHEzAGdGY3nkit7bwTZUO7NcffUwy5WrYrmOftmOUMgun5OvUM8LVzmbqvY5WeFqfeHTc+SLnOEKI2bXMJ2mpLSE0m8LKycv18nrN7IdNvQVUgCWUI7RClgKwBFSAI5Q1dIxpAAcIeniGlIArtC3lBSANSy812HDF0Q4VNhCoTs7a6g4KnwHbMP3w8rrr6xXjn6FbcNxT7n6N6QMVf10V/vnkXu/f/n/us6rkZKm8c+saxBV3AAgk8kAAH6NF4u/QJ92Mbe/62mViIjIg0wmg0QiUZZlR6NRpFIp/Dqt+V0xCalUCtHo+PLzf+0sVWEnQFzXxVtvvYW6ujpYloW+vj7Mnj0bXV1dqK+vP9erV7G4nSaG22liuJ0mhtupOKUUMpkMZs2aBdsu35nm4eFhZLPyLMVERKNRxOP6e738taq4GQDbtnHRRReNe76+vp5fsAngdpoYbqeJ4XaaGG6n8cr1l/8HxePxQP7iniq8CJCIiCiAOAAgIiIKoIofAMRiMTz88MOIxWLnelUqGrfTxHA7TQy308RwO5GfVdxFgERERFR+FT8DQERERFOPAwAiIqIA4gCAiIgogDgAICIiCqCKHwA89thjaGhoQDwex8KFC/HKK6+c61U6p15++WV84QtfwKxZs2BZFv7zP/9zTLtSCs3NzZg1axaqqqqwbNkyvPHGG+dmZc+RlpYWfOITn0BdXR1mzJiBW2+9FceOHRvzGm4nYPv27bjqqqsKN7FZvHgxfvnLXxbauY2Ka2lpgWVZWLNmTeE5bivyo4oeADz33HNYs2YN1q9fj8OHD+P666/HTTfdhDfffPNcr9o5MzAwgKuvvhqtra1F2zdv3owtW7agtbUVBw4cQCqVwvLlyws1FoKgo6MD99xzD/bv34/29nbk83k0NjZiYGCg8BpuJ+Ciiy7CI488goMHD+LgwYP41Kc+hVtuuaXwi4vbaLwDBw5gx44duOqqq8Y8z21FvqQq2LXXXqvuvvvuMc9ddtll6rvf/e45WqPKAkDt3r278H/XdVUqlVKPPPJI4bnh4WGVSCTU448/fg7WsDJ0d3crAKqjo0Mpxe0kOf/889WPf/xjbqMiMpmMmjt3rmpvb1dLly5VDzzwgFKK+xP5V8XOAGSzWRw6dAiNjY1jnm9sbMTevXvP0VpVtuPHjyOdTo/ZZrFYDEuXLg30Nuvt7QUATJs2DQC3UzGO42DXrl0YGBjA4sWLuY2KuOeee/C5z30On/nMZ8Y8z21FflVxxYDOOHXqFBzHQTKZHPN8MplEOp0+R2tV2c5sl2Lb7MSJE+dilc45pRTWrl2LT37yk5g/fz4AbqcPOnLkCBYvXozh4WHU1tZi9+7duOKKKwq/uLiNRu3atQuvvfYaDhw4MK6N+xP5VcUOAM6wLGvM/5VS456jsbjN/uLee+/F66+/jl//+tfj2ridgHnz5qGzsxOnT5/G888/jzvuuAMdHR2Fdm4joKurCw888AD27NkjVp7jtiK/qdhTANOnT0coFBr31353d/e4kTaNSqVSAMBt9r777rsPL7zwAn71q1+NKTHN7fQX0WgUH/3oR3HNNdegpaUFV199NX74wx9yG33AoUOH0N3djYULFyIcDiMcDqOjowM/+tGPEA6HC9uD24r8pmIHANFoFAsXLkR7e/uY59vb27FkyZJztFaVraGhAalUasw2y2az6OjoCNQ2U0rh3nvvxc9//nP8z//8DxoaGsa0czvpKaUwMjLCbfQBn/70p3HkyBF0dnYWHtdccw2+9rWvobOzE5dccgm3FfnTubv+0GzXrl0qEomoJ598Uv3+979Xa9asUTU1NepPf/rTuV61cyaTyajDhw+rw4cPKwBqy5Yt6vDhw+rEiRNKKaUeeeQRlUgk1M9//nN15MgR9dWvflXNnDlT9fX1neM1P3u+9a1vqUQioV566SV18uTJwmNwcLDwGm4npdatW6defvlldfz4cfX666+rhx56SNm2rfbs2aOU4jaSfDAFoBS3FflTRQ8AlFJq27Zt6uKLL1bRaFR9/OMfL0S5gupXv/qVAjDucccddyilRiNJDz/8sEqlUioWi6kbbrhBHTly5Nyu9FlWbPsAUE899VThNdxOSn39618vfLcuvPBC9elPf7rwy18pbiPJhwcA3FbkRywHTEREFEAVew0AERERlQ8HAERERAHEAQAREVEAcQBAREQUQBwAEBERBRAHAERERAHEAQAREVEAcQBAREQUQBwAEBERBRAHAERERAHEAQAREVEAcQBAREQUQP8fM6BE1aZEUxoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(mass_collision)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
